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Background: Over the last few years novel folding mechanisms of globular proteins have been 
proposed using minimal lattice and off-lattice models. However, the factors determining the coop- 
erativity of folding in these models and especially their explicit relation to experiments have not 
been fully established. 

Results: We consider equilibrium folding transitions in lattice models with and without side 
chains. A dimensionless measure, f2 c , is introduced to quantitatively assess the degree of coopera- 
tivity in lattice models and in real proteins. We show that larger values of f2 c resembling that seen 
in proteins are obtained in lattice models with side chains (LMSC). The enhanced cooperativity 
in LMSC is due to the possibility of denser packing of side chains in the interior of the model 
polypeptide chain. We also establish that O c correlates extremely well with a = (Tq — Tf)/Tg, 
where Tq and Tf are collapse and folding transition temperatures, respectively. These theoretical 
ideas are used to analyze folding transitions in various two state folders (RNase A, chymotrypsin 
inhibitor 2, fibronectin type III modules, and tendamistat) and three state folders (apomyoglobin 
and lysozyme). The values of £l c extracted from experiments show a correlation with a (suitably 
generalized when folding is induced by denaturants or acid). For lattice models and for proteins 
Sl c ~ exp[-(o-/o- ) /3 ], where /3 < 1. 

Conclusions: A quantitative description of cooperative transition of real proteins can be made 
by lattice models with side chains. The degree of cooperativity in minimal models and real proteins 
can be expressed in terms of the single parameter a, which can be estimated from experimental 
data. 

Key words: Lattice models with side chains, cooperativity of folding, collapse and folding 
transition temperatures, chemical and thermal unfolding. 

I. INTRODUCTION 

In recent years, considerable insight into thermodynamics and kinetics of protein folding 
has been gained using simple lattice and off-lattice models [0-0. The simplicity of these 
models allows for exhaustive sampling of the available conformational space, so that a de- 
tailed picture of the folding mechanisms can be obtained. In almost all these studies only 
a very coarse grained description of the polypeptide chain is retained. These models can 
be viewed as providing, in an approximate way, a simple representation of the a-carbons of 
proteins. The amino acid residues are represented as single beads, which is very much in 
spirit of united atom approach. Mapping the characteristics of such a coarse grained descrip- 
tion onto real polypeptide chains suggests that each bead may correspond to 2-3 amino acid 
residues |7|,|26|| . Lattice models (LM) further restrict positions of the beads to the vertices of 
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an appropriate lattice. With the exception of the works by Skolnick and coworkers || most 
of the lattice studies have utilized simple versions of cubic and square lattices. 

Despite the success of the minimal models of proteins in providing general principles 
of folding mechanisms there are several issues that have to be resolved before a realistic 
description of folding can be provided. In this article, we confine ourselves to one issue, 
namely, to what extent is the cooperativity of folding from an unfolded state to the native 
conformation mimicked by various minimal models. It is known experimentally that un- 
der optimal folding conditions proteins (especially small ones) fold in a highly cooperative 
manner, i.e., they exhibit an all-or-none transition ||. This implies that there are only 
two relevant thermodynamic states, namely, the unfolded state (U) and the native folded 
conformation (N). The transition U # N is a first order like phase transition, implying that 
there is a barrier (often modest) separating the two states. The first order nature of folding 
transition has been seen in various minimal models of proteins for a variety of interaction 
schemes. Furthermore, it has been shown by various groups that this transition is a coop- 
erative |l|-|6|,p!0|-p^| depending on a number of factors, which are intrinsic to the sequence. 
However, it has not been established the degree to which this cooperativity mimics that seen 
in real proteins. 

To address this issue quantitatively, we study lattice models with and without side chains 
with the purpose of assessing the extent to which cooperativity of folding transition in real 
protein is mimicked by these models. We introduce a dimensionless measure of cooperativ- 
ity that allows us to compare theoretical models and data on real proteins on nearly equal 
footing. Previous studies have also addressed the factors that determine the cooperativity of 
folding using simplified models [§,[L4|,|nj . Hao and Scheraga showed, using a model similar to 
that employed here, that cooperativity of folding transitions results from proper long ranged 
interactions |T3| , |I5|] . More recently, Skolnick and coworkers have used a high coordination 
lattice with each side chain having many rotamer degrees of freedom ||. They showed that 
cooperativity, as measured by the sharpness of the specific heat curves, results from tight 
side chain packing [pl|[17]]. These results together with our analysis using the cooperativity 
measure clearly show that to achieve realistic description of thermodynamics one has to 
include the rotamer degrees of freedom (at least approximately) into coarse grained mod- 
els so that dense core packing seen in real proteins can be simulated. We establish that 
cooperativity can be linked to experimentally measurable quantities, namely, the collapse 
transition temperature Tq and the folding transition temperature Tf. Both Tq and Tf are 
dependent on sequence, pH, and concentration of denaturants. 



II. METHODS 

A. Lattice Model with Side Chains (LMSC) 

A protein with side chains is modeled on a cubic lattice by a backbone sequence of 
iV beads, to which a "side" bead, representing a side chain, is attached. This method of 
incorporating the rotamer degrees of freedom in lattice models is similar to that proposed 
by Bromberg and Dill |18 |. In this coarse grained description, the peptide bond and the 



a— carbon are given by a single bead. The system has in total 2N beads. The following types 
of interactions are included in this model, (i) Self-avoidance, i.e. any backbone and side 
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beads cannot occupy the same lattice site more than once, (ii) The interactions between side 
beads % and j — > 1) are taken to be pairwise and are assumed to be given by statistical 
potentials derived by Kolinski, Godzik, and Skolnick (KGS) (ref. ||19|| , Table III). We assume 
that, if two side beads i and j (\j — i\ > 1) are nearest neighbors on a lattice, they form 
a contact and the energy of their interaction is given by the KGS potential matrix element 
associated with the residues i and j. The interaction energies in the KGS interaction matrix 
are in units of kT, where T is the absolute temperature. In two dimensional representation 
we have also tried Miyazawa-Jernigan (MJ) |20| , |2T| parameterization of the matrix of contact 
interactions between amino acids as well as random site interaction model. 

In this study, we took N = 15. It is known that the total number of symmetry unre- 
lated self-avoiding conformations on a cubic lattice (i.e., for a model without side chains) is 



93,250,730 for iV = 15 [P2" j. Introduction of side beads dramatically increases the available 
conformational space, and hence exact enumeration of all conformations is no longer feasible 
even for iV = 15. To illustrate the sharp increase in the number of available conformations 
we take square lattice as an example. On a square lattice the total number of symmetry un- 
related conformations is 296,806 for N = 15 fllQ |. When side beads are attached the number 
of available conformations becomes 653,568,850 and the corresponding number of backbone 
conformations is 204,988. The number of backbone conformations is reduced because many 
compact backbone conformations are now inaccessible due to side chain overlaps. There- 
fore, inclusion of side chain rotamers in the model increases the number of conformations 
by a factor of 2202. In the case of cubic lattice this factor would be even larger because 
the coordination number for cubic lattice is 6 instead of 4 for a square lattice. Thus, one 
has to resort to Monte Carlo methods to obtain thermodynamic characteristics of sequences 
even for N = 15 for cubic lattice models with side chains. The maximum number of near- 
est neighbor contacts c in a fully compact structures also increases, when side chains are 
included. In our our previous work |22] we showed that for N = 15 c is 11 in the model 
without side chains. For models with side chains included this number increases to 22, when 
only the contacts between side chains are counted. The maximum total number of nearest 
neighbor contacts between backbone beads as well as side chains is 30. The larger value of c 
gives rise to a denser packing of the interior, and hence to a greater degree of cooperativity. 



B. Selection of Sequences and Determination of Native Structures 

We generated, at random, several sequences using the constraint that the fraction of 
hydrophobic residues in a sequence be approximately 0.5 as in the case of naturally occurring 
proteins. Random sequences are most likely to be not protein-like. In order to obtain viable 



protein-like sequences we utilized the Z score optimization of random sequences [[23| -|25 
having unique ground states. The idea of Z score optimization is based on minimizing 
the relative difference between the energy of a sequence in the target conformation E Q and 
the average energy of misfolded structures E ms , an estimate for which is E ms = c< B >, 
where c is number of nearest neighbor contacts in misfolded structures (presumably equal 
to the number of contacts in the native conformation) and < B > is the average contact 
energy calculated for a particular sequence. Thus, the quantity minimized in the course of 
optimization procedure is 
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z- - s (1) 

where 5 is the standard deviation of contact energies for a given sequence. Standard deviation 
S and E ms are sequence dependent quantities. Low values of Z would correspond to sequences 
with large stability gaps (E ms — E ) and small 5. The sequence with lowest Z score for 
a given sequence composition is found by performing Monte Carlo simulations in sequence 



space, the details of which are given elsewhere p2p4"p5f . In this study, we used simulated 
annealing protocol when doing Monte Carlo simulations in sequence space and we produced 
a total of 20 trajectories. In the course of simulations we constantly monitored Z score 
values and in the end selected sequences with low values of Z. 

Target conformations used in optimization procedure were randomly chosen from prelim- 
inary Monte Carlo simulations subject to the conditions that these structures be reasonably 
compact, and have a unique contact map (in order to avoid degeneracy of a ground state). 
When a sequence was optimized, we verified that the target conformation indeed corre- 
sponds to the lowest energy state. To this end, we generated 20 slow cooling Monte Carlo 
trajectories starting with different initial conditions at high enough temperature, at which 
sequences are in random coil conformations. The final temperature of these simulations is 
set to 0.05, which is well below the folding transition temperature, Tf. We were unable to 
find any conformation, which has lower energy for selected sequence than the target struc- 
ture. We also constantly monitored energies in all the dynamic simulations. In no instance 
did we find energies lower than that corresponding to the native state. The energies and 
temperature in this study are expressed in the units of kT and T, respectively, where T is 
an absolute temperature. 



C. Probes of Thermodynamics of Sequences: Computation of Tg and Tf 

Thermodynamics of a given sequence was probed using the overlap function \ 
and its fluctuation A% suitably modified for the lattice model with side chains. The overlap 
function, which describes the overall conformation of a sequence in relation to the native 
state N is defined as 

y = 1 - 2^-3JV + l (^^- r g ,J,)+ E't'S-'D+IW "'if')). < 2 > 
ziv oiv -t- i iKj i<j+1 ^ 

where r|j , rfj, rfj refer to the distances between backbone beads, between side beads, and 
between backbone and side beads, respectively. The factor 2N 2 — 3iV + 1 ensures that in 
the native conformation \ — 0- We also studied the partitioned overlap functions, Xss and 
Xbb, which are defined by 

X« = 1-]^E^-^) (3) 

i<j 



and 

Vll =1 

N 2 — 3N + 2 



X»=l- vrL £ S($-r$"). (4) 

i<j+l 



4 



The backbone overlap Xbb measures the content of a native backbone, while the side chain 
overlap Xss measures the "nativeness" of side chain conformations. By treating these quan- 
tities separately it is easier to monitor the onset of native character in the backbone and 
side chains. 

The acquisition of the overall native state (backbone plus side chains) was also probed 
using the probability of being in the native structure (a single microstate) 



El 

T 



Po = E ^ (Xv ; ' • (5) 

In order to observe the acquisition of the native conformation of the backbone and side beads 
separately, we have calculated the probability that the backbone is in the native structure 

p* = ^5( Xhbb )e-^ 
and, similarly, the probability of side chains to be in the native positions 

° ~~ ~ ns • ' ' > 

z2i e T 

It may seem that there is no need to consider P S separately from P , because there are not 
many ways of changing the backbone conformation while keeping the side chain contacts 
intact. Indeed, analysis shows that there are very few conformations, which do allow such 
changes, but they constitute a small fraction of all structures so their contribution would be 
negligible. Our simulations show that unless these structures have a considerable Boltzmann 
weight (that is highly unlikely) there is no appreciable difference between Po and P S . 

In Eqs. (|§,§,0) the various probabilities correspond to one microstate associated with 
the native conformation, which is an artifact of lattice models. This is physically unrealistic 
because in any realistic description there is a volume associated with the native basin of 
attraction (NBA). Thus, many conformations will map onto the NBA. It is, therefore, more 
relevant to calculate the probability of being in the NBA 

r^NBA — — gr • (8 ) 

In Eq. (|8[) the sum is taken over all possible conformations of backbone and side chains, and 
Xnba is the value of overlap at the temperature, when the fluctuations of x are maximum, 
i.e. Xnba =< x(Tf) > (see below). At and below Tf conformations with x < Xnba "will map 
onto the NBA. It is clear that Pnba is physically more relevant than P . This is because 
thermal fluctuations in real proteins frequently sample numerous conformations proximal to 
the native conformation. 

In addition to the functions described above, we calculated the energy <E>, the specific 
heat C v , and the radius of gyration R g using standard definitions. The collapse transition 
temperature Tg was determined from the maximum of specific heat C v HlO|, 22 ,p6|| and was 



found to be well correlated with the temperature, at which the rate of change of the radius 
of gyration R g becomes maximum. As before the folding transition temperature Tf was 
inferred from the peak of the overlap fluctuation Ax [0,^,^] . 
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D. Calculation of thermodynamic quantities 



The thermodynamic quantities were computed using the multiple histogram technique 
27], which has been implemented in the following way. From the slow cooling simulations 
used to determine the native conformations a rough estimate of collapse and folding transi- 
tion temperatures, Tq and Tf, can be made. This allows us to select the temperature range, 
in which histograms are to be collected. At each temperature histograms were collected 
using M = 100 — 200 Monte Carlo trajectories, each with different initial conditions. Each 
trajectory starts at high temperature T h > T e and ends at the temperature TJ ~ Tf. In 
the course of a trajectory the temperature is changed every t Monte Carlo steps (MCS) 
according to the schedule T r = Th — rAT, where r = 0, {Th — T\)j AT and the typical 
value of AT is 0.05. Accordingly, histograms were collected at the temperatures T r using M 
trajectories generated. In total, we have R = (Th — Tj)f AT + 1 histograms. Five variables 
were used to collect histograms, these are overlap function x, backbone overlap Xbb, side 
chain overlap Xss, energy E, and the square of the radius of gyration R 2 g . The length of the 
trajectory at fixed temperature T r was set constant for all temperatures and sequences at the 
value of t = 2 x 10 7 MCS. It is important to note that a portion of trajectory immediately 
following temperature change must be excluded in order to allow the chain to equilibrate at 
a new T. We found this interval of equilibration does not exceed 2 x 10 6 MCS at the lowest 
temperature for all sequences. For simplicity, the equilibration interval was kept constant 
for all temperatures. While collecting states in the histograms we used 0.1 grid intervals 
for the energy E and the square of the radius of gyration R 2 and exact discrete values for 
functions x, Xbb, and Xss- 

The thermodynamic average of a quantity A that is a function of, say, x an d E is 
calculated using multiple histogram method as follows 



Er=l ™" e r 

where R is the number of histograms, T r is the temperature of simulations, at which rth 
histogram is collected, n r is the number of states in the rth histogram, and f r is the scaled 
free energy to be calculated self-consistently from the following equation 



' = vv e -l^!M^ (io) 



The values of f r can be determined from an iterative scheme within less than 300 iterations 
with excellent accuracy. 

Using the multiple histogram technique we have calculated overlap function <x>, its 
fluctuation A%, backbone overlap <Xbb>, its fluctuation Axbb, side chain overlap <x ss >, 
its fluctuation Ax ss , energy <E>, specific heat C v , probabilities of occupancy of the native 
state Po, Pq, Pq and the NBA Pnba, and the radius of gyration < R g > as functions of 
temperature. From these the key characteristic temperatures Tq and Tf can be readily 
computed. 



6 



III. RESULTS AND DISCUSSION 



A. Thermodynamics of sequence A and B 

For the present study we have considered three sequences in three dimensional (3D) 
LMSC each with a unique ground state. In order to keep the description compact we 
describe in detail the thermodynamics associated with two sequences, A and B. The general 
characteristics of the third sequence C are similar to those of B. The compositions of the 
two sequences along with their native states and energy spectra are given in Fig. (1). 

In Fig. (2) temperature dependencies of various thermodynamic functions are plotted 
for sequence A. The collapse transition temperature T 9 is determined from the peak of C v 
(see Fig, (2c)) and is equal to 0.27. As shown in Fig. (2c) this temperature coincides 
with the temperature, at which the derivative d < R g > jdT shows maximum. This is a 
strong evidence that T e indeed corresponds to the temperature, at which the transition 
from a random coil state to an ensemble of compact structures takes place. The plot of 
< x > (measuring the content of native state) shows that equilibrium refolding occurs highly 
cooperatively, indicating a two state transition. Note that the plot of < Xbb > follows <x> 
very closely. It can be seen in Fig. (2b) that the fluctuations in the functions x an d Xbb (as 
well as Xss, not shown) display a well defined peak at the temperature Tf = 0.26. This value 
coincides with the temperature, at which the probability of being in the NBA is 0.5 (see 
Fig. (2d)). Thus, evaluations of Tf from Ax and Pnba are equivalent. It is seen that the 
backbone fluctuations Axbb are stronger at Tf than the fluctuations of the entire sequence. 
Thus, for sequence A the value of parameter ot is 0.04 (see below). 

It is of interest to compare the probabilities Pq and Pq as a function of temperature. 
We find that the native backbone is formed prior to the formation of the entire native 
conformation as the temperature is lowered (data not shown). In fact, the native backbone 
is predominantly populated at T < T b = 0.22, whereas the complete native structure 
becomes overwhelmingly populated only at T < 0.17. This observation suggests that as the 
temperature decreases backbone conformation becomes frozen, while side chains still have 
no well defined conformations seen in the native structure. 

Our analysis indicates that sequence A is a two-state folder. This is readily inferred from 
Figs. (2a) and (2b). This point is further illustrated by plotting the thermally weighted 
distribution of (x, E) states at T ~ Tf. Such plots can be readily obtained from multiple 
histogram calculations. In Fig. (3) we plot the two dimensional histogram of (x, E) states 
for sequence A (upper panel). The main feature of this plot is that it clearly reflects the 
localization of sequence in two distinctive states shown in red with very different values of 
overlaps and energies - folded (NBA) and unfolded states. Kinetic simulations (D.K. and 
D.T., unpublished) suggest that the second (unfolded) state consists of open random coil 
conformations. Therefore, we classify sequence A as a two state folder. 

The thermodynamic analysis for sequence B gives T e = 0.44. At roughly the same 
temperature (=0.48) d < R g > / dT has a maximum indicating that Tq does correspond to 
the collapse transition. The value of Tf = 0.30 and hence a = 0.32. Due to this relatively 
large o value we expect the folding transition to be better described by a multistate behavior. 
This is in fact the case and can be seen most clearly from the two dimensional histogram of 
(x, E) states displayed in the lower panel of Fig. (3). This plot for sequence B shows that, 
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at the folding transition temperature, the model polypeptide chain samples at least three 
distinctive thermodynamic states shown in red with x — 0.0, x — 0.3, and x — 0.65. Notice 
that the energies at these values of x are a l so different. Based on this analysis we expect 
the folding of B to be less cooperative than A. 



B. Comparing cooperativity for LMSC and LM 

The structures of the native states and the corresponding low energy spectra for the two 
sequences A and B are given in Fig. (1). Even though the chain length is small it is seen 
that an interior of the native state of a sequence A is more tightly packed than that of B. By 
contrast, sequence B has less well formed interior and hence has a lesser degree of side chain 
packing. Specifically, the radius of gyration of a hydrophobic core R 9j h-h (which measures 
the average distance between hydrophobic residues) calculated for the native structure of 
sequence A is equal to 1.2. In contrast, R 9 ,h-h obtained for sequence B is found to be 1.25. 
It is also interesting to point out that the average distance between hydrophilic residues in 
the native structure of sequence A is considerably larger (equal to 2.3) when compared with 
that found for sequence B (equal to 1.9). This suggests that the core of sequence A has 
significant hydrophobic clustering that leads to denser packing than in sequence B. The low 
energy spectra of these sequences also reveal different behavior. The first excited state of 
sequence A has the backbone in a non-native conformation, while the first few excited states 
of B all have the native conformation for the backbone. The energy gap (the difference 
between the native state and the first excited state) and the so called stability gap |2J (the 
energy difference between the native conformation and misfolded structures) are larger for 
sequence B than for sequence A. Nevertheless, it is shown below that sequence A folds more 
cooperatively than sequence B indicating that the gaps (however they are defined) are by 
themselves not a good indicator of cooperativity. 

In order to compare the degree of cooperativity in various models and real proteins we 
introduce the dimensionless cooperativity index defined as 

max\ d<x> ] 

"c — J-max ^rp i \ LL ) 

where AT is the full width at the halfmaximum of the peak of d < x > /dT and T max is 
the temperature, at which d < x > /dT has a peak. This quantity can be calculated for 
other measures such as the temperature or denaturant dependence of the fraction of native 
state, which are experimentally measurable. 

There are other measures of cooperativity that could be used. One such criterion is 
based on the van't Hoff equation, from which an index of cooperativity, namely, &vh — 
AHAT/RT^ ax may be introduced. For an infinitely sharp two-state transition ttvH — * 
const, while for a completely non-cooperative transition Qvh — * 0. Our criterion, on the 
other hand, gives Q c —>■ oo for an infinitely sharp two-state transition and tends to zero, 
when the transition is non-cooperative. This provides a wider spread in Q c and we expect 
it to lead to easier calibration of cooperativity. Furthermore, when only the dependence of 
order parameter (like fraction of native state) as a function of denaturant concentration is 
available, a suitable generalization of Eq. ([□]) can be readily used to compute Q c , whereas 
the generalization of flyH to this case is not obvious. 
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We have shown elsewhere [TO,22|,26,p8| that folding kinetics (i.e., the time required to 



reach the native conformation) is well correlated with 

t = — 7f , (12) 

i.e., sequences with small a reach the native state more rapidly than those with larger a. 
Folding is most rapid near the tricritical point (Tf w Tg), at which collapse and acquisition 
of the native state are almost synchronous. Parenthetically, we note that other criteria for 
identifying fast folding sequences have been proposed f2|j6[]. 

The value of Q c for sequence A is 5.32, for sequence B it is 2.03, and for the third 
sequence, C, Q c = 1.80. Comparing Q c for sequences A, B, and C shows that cooperativity 
is also determined by a; a for A is 0.04, for B is 0.32, and for C is 0.34. This is shown (see 
below) to be the case for both LMSC and LM. 

It is interesting to compare Q c for lattice models with side chains and lattice models 
without side chains. For this, we consider a number of 15-mer sequences without side chains 
on a cubic lattice. The sequences used for this purpose are the ones studied in our previous 
papers p2[28|. Many of these sequences fold cooperatively with a values ranging from 0.0 to 



0.79. In Fig. (4) we plot Q c as a function of a for 15-mer lattice sequences. The maximum 
value of Q c of 3.0 is observed for sequences with a ~ 0.0. This level of cooperativity is less 
than that seen in the models with side chains indicating that cooperativity in natural proteins 
results from the dense packing of side chains in the core |8]JT6|]. Fig. (4) shows remarkable 



statistical correlation between the cooperativity index f2 c and a. When combined with our 
previous results we can conclude that both cooperativity and folding kinetics is determined 
by the intrinsic thermodynamic parameters Tg and Tf. It should be emphasized that Tg and 
Tf are sensitive to mutations as well as external conditions such as pH, ionic strength, etc. 

In order to further establish the correlation between Q c and a we have studied two di- 
mensional (2D) LMSC using a number of sequences with very different interaction potentials 
(see caption to Fig. (5)). Sequence properties have been inferred by exhaustive enumeration 
of all allowed conformations. This permits us to analyze much larger sequence database 
than for 3D LMSC. In Fig. (5) we plot Q c s a function of a for our 2D models with side 
chains. As is the case for N = 15 3D LM the correlation is excellent. It should be noted 
that Q c values are considerably smaller than in 3D indicating that for realistic depiction of 
cooperativity 2D simulations may not be reliable. 



C. Dependence of £l c on chain length 

The model with side chains for iV = 15 yields Q c values of the order of 5, when a is small, 
whereas real proteins, which are two state folders, have Q c in excess of 10 (see below). This 
discrepancy between LMSC and experiments is due probably to the small value of N. If N is 
increased, then we expect to obtain Q c in the experimentally measured range. This assertion 
is borne out by examining Q c for lattice models without side chains for larger N. In Fig. 
(6) we display fi c vs a for 27-mer sequences, for which detailed account of thermodynamics 



and kinetics was provided in a previous study For iV = 27 we see that Q c ranges from 
1.0 (a = 0.35) to 33.3 (a = 0.0). This figure also shows that cooperativity is essentially 
determined by a. The results in Figs. (4-6) show that one can write 
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tt c ~ exp[-((x/a )% 



(13) 



where o"o is a model dependent constant and j3 ^$ 1. Since the addition of side chains enhance 
fl c , we conclude that in order to provide a realistic description of cooperativity in protein 
folding it is necessary to include side chains in minimal models of proteins. 



D. Extracting Sl c and a from experiments 

We have computed fl c for a few proteins and a protein fragment using experimental data 
such as the dependence of fraction of native state as a function of temperature, denaturant or 
acid [p9| -|35H . The proteins we chose for our analysis are for illustrative purposes only and does 



not constitute an exhaustive analysis. The data from the literature when used in combination 
with Eq. (|TT| ) gives Q c in the range 0.25 - 99.4 (see Tables I-III). Out of the eight proteins 
considered here the fragment of barnase is least cooperative []3"3| . This, in fact, suggests that 
the cooperativity is greatly compromised in the absence of long range interactions with the 
rest of the protein. Surprisingly, the cooperativity of the ninth module of fibronectin exhibits 
the level of cooperativity similar to the barnase fragment. Chymotrypsin inhibitor 2 (CI2), 
RNase A, tendamistat, and the tenth module of fibronectin ( 10 FNIII) have been recently 
shown |29-p2| to fold in a very cooperative two state manner. Apomyoglobin (apoMb) 



and lysozyme |$4],|35l have been shown to be three state folders. Below we describe the 



cooperativity of these proteins using the equilibrium folding data. 

Two state folders: The thermal unfolding of RNase A and RNase B has been recently 
reported |3*0|| . The data for the fraction of native state as a function of temperature (see 
Fig. (5) of [30p) clearly exhibits a cooperative two state transition. Using the parameters 



obtained from the two state analysis of RNase A in [|3D| and Eq. ([□]) (where < x > 1S 
replaced by the native fraction f^) we get Q c ~ 99.4. The large value of fl c makes us 
suggest that for this protein a ~ 0.0 implying that Tf « Tq. A more precise estimate of o 
is given below. 

Kiefhaber and coworkers [31] have argued using the ellipticity measurements at pH=7 



and T = 25°C that tendamistat (a small protein consisting of 74 amino acid residues) 
unfolds thermodynamically in an all-or-none fashion upon increasing the concentration of 
guanidinium chloride. The value of fl c calculated using the dependence of the fraction of 
native state fa (see Eq.(|i~4"D) on the concentration of GdnHCl is 13.7, which is considerably 
less than for RNase A. This suggests that either a for tendamistat is larger than for RNase 
A or that fl c is intrinsically smaller when folding is induced by varying the concentration of 
denaturants as opposed to thermal unfolding. The determination of a (see below) suggests 
that the latter is a more likely possibility. 

Recently, Dobson and coworkers have shown that denaturation induced unfolding of 
fibronectin type III 9 FNIII and 10 FNIII proteins is well described by a two state model 
32| . However, it was found that 9 FNIII has a lower stability (as measured by a free energy 



difference between U and N) than 10 FNIII. The free energy of unfolding in water AGh 2 o — 
1.2 and 6.1 kcalmol~ l for 9 FNIII and 10 FNIII, respectively. In line with these observations, 
our analysis suggests that both proteins display sharply different degrees of cooperativity 
with Q c = 0.28 obtained for 9 FNIII and 7.5 for 10 FNIII. Accordingly, we expect that the a 
values are very different for 9 FNIII and 10 FNIII. 
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A few years ago Jackson and First p9j reported that folding of CI2 proceeds via a two 



state manner. This was the case whether denaturation was induced chemically or thermally. 
By using fluorescence, F, as a function of [GdnHCl] they extracted parameters for the two- 
state description of the data. Since the values of F^ and Fy (fluorescence in the native state 
and denaturated state, respectively) were not reported we determine the fraction of native 
state for chemically induced unfolding using 

fa = ~p V~ = ~ SoF' ( 14 ^ 

fu-fN l + exp"-ra^ 

where AGu = AG H 2 o— m [GdnHCl] and AGh 2 o is the free energy of unfolding in water, m is 
a constant that effectively measures the degree of exposure of the protein upon denaturation, 
and R is the gas constant. For Eq. ([14]), which gives a two-state description of denaturant 
induced unfolding, fl c can be shown to be 




q c = ii — - (15) 

" i 3-2v / 2 

implying that cooperativity in this case is directly related to stability fl3T |. Using the ex- 



perimental parameters AGh 2 o — 7.03 kcalmoF 1 and m = 1.79 kcal mol~ l M~ x |2pJ we find 
Q c = 9.9 - a value that is comparable to tendamistat. 

Jackson and Fersht also reported thermally induced unfolding of CI2 as a function of pH 
29fl . They measured the pH dependence of the folding transition temperature Tf (referred 



to as T m in ref. ||29|| ) using scanning microcalorimetry experiments. They showed that in the 
range of pH (2.2-3.5) CI2 unfolds thermally also by a two state manner. We calculated Q 



c 



using the data given in [2D] and find that it lies between 5.9 at pH=2.2 to 29.3 at pH=3.5. 
The degree of cooperativity clearly increases with pH. At the higher values of pH Q c is larger 
than that found by chemical denaturation. Thus, in general, it appears that the degree of 
cooperativity is greater when unfolding is induced thermally than by adding denaturants. 

We have established in the simulations on lattice models that fl c is well correlated with 
<jt (see Figs. (4-6)). In order to verify this in experiments one needs measurements of Tg, 
which are not currently available. So we have devised a way to obtain Tg using the following 
arguments. Since there are two separate phase transitions at Tq and Tf the experimen- 
tal determination of a requires separate measurements of order parameters, probing these 
"phases". The estimate of T g requires independent measurement of a quantity associated 
with the compaction of a polypeptide chain (for example, the radius of gyration, R g ). At 
present, simultaneous measurements of R g and the fraction of native state Jn are extremely 



rare |35|. In the absence of such measurements we have used an ad hoc empirical procedure 



to estimate a from the existing data so that the variation of Q c with pH, such as that seen 
in CI2, can be rationalized. To estimate Tg from experiments we propose the following 
approach. We calculated Tq (an estimate for Tg) using the condition 

fam = c, (i6) 

where C here is taken to be 0.05. At Tg the equilibrium is shifted almost to the unfolded 
state and, consequently, we expect C in Eq. (|1^) to be small. Recall that Tf is determined 
in the analysis of experimental data from f^(Tf) = 0.5. 
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In order to justify the use of Eq. (0) we have estimated T' e using Pq(Tq) = 0.05 for LM 
with iV = 27. A plot of T' d from this estimate as a function of actual Tq (from the peak 
of the specific heat) is given in Fig. (7). There is a very good correlation (with the slope 
being near unity) with the quality getting better as Tq becomes larger, which corresponds 
to small cr values. A similar correlation is found when C is changed to 0.1. This method 
of estimating Tq, although arbitrary, is physically motivated. We find that if C is increased 
beyond 0.1 the correlation coefficient considerably decreases. Thus, this way of obtaining 
Tq from experimental data (with C ;$ 0.1) is useful especially for the two state folders. 

Additional justification of the use of Eq. (|T6|) to estimate Tq can be given by examining 
the temperature dependence of R g for RNase A. Sosnick and Trewhella have used small 
X-ray scattering to measure R g as a function of temperature (see Fig. (2) of fl36H ). For 
the nonreducing conditions, which were the ones considered in [|30| as well, we estimate Tq 



(from Fig. (2) in J36|) to be around 66°C. This value is remarkably close to that obtained 
using Eq. flT6"D (Tq ~ 65.2°C) using the data of Ulbrich-Hofmann and coworkers. We 
should emphasize that despite the reasonable estimates of Tq from Eq. ([16]) independent 
experimental measurements are needed. 

We have used Eq. (|16"D together with the data on thermal unfolding of CI2 and RNase A 
to estimate Tq. These are given in Table I along with their Q c and a. It is seen from Table 
I that ot is relatively large for CI2 at pH=2.2. From the previously established results 



connecting folding kinetics and a p2| , |28| , |26|1 we predict that there should be a significant 
off-pathway process in the refolding of CI2 at pH=2.2 at T pa Tf. It would be interesting 
to determine the pH dependence (over a wider range) of Tf for the variant of CI2 used in 
the recent experiments of Fersht and coworkers, which contains 64 amino acid residues as 
opposed to 81 in the original work of Jackson and Fersht. 

In most of the experiments the equilibrium unfolding is achieved by chemical denatura- 
tion. The analysis of such experiments can be done by generalizing a. We introduce 

°d = — ^ , (17) 

where Cq is the concentration of denaturant, at which the transition to compact structures 
from the random coil state takes place, and Cf is the denaturant concentration in the 
midpoint of an appropriate measure of the native state population. Similarly, Q c in Eq. 
(PI) can be calculated using the denaturant dependence of the fraction of native state. 
With this generalization we have estimated Cf and Cq for tendamistat, 9 FNIII, 10 FNIII, 
and CI2 (denaturant induced unfolding). The value of Cq is obtained from the condition 
In(Cq) = 0.05. The values of cr and Q c are given in Table II. 

The physical basis for proposing Eq. ( pj) as a natural generalization of Eq. ( |T2"D is 
the following. When the temperature changes from below Tf to above Tq, there are two 
phase changes, one at Tf and the other at Tq. Similar phase changes are seen upon chemical 
denaturation. Both a D and cr^ capture the physics that fastest folding should occur when the 
intermediate is eliminated. Quantitatively, the value of <jd and <jt need not coincide because 
the nature of interactions due to thermal unfolding and denaturant induced unfolding is 
different. In fact, it appears from the analysis of the experiments in this study that thermal 
denaturation is sharper than chemically induced unfolding. 

Three State Folders: Barrick and Baldwin [R3] showed that urea and acid induced folding 
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of apomyoglobin (apoMb) is best described by a three state model. Thus, the equilibrium 
folding is given by 



U^I^N. (18) 

Here we analyze the acid induced transition using the data presented in Fig. (1) and Table 
I of [ 3"3fl . In Fig. (1) of Barrick and Baldwin present the helical content of apoMb 
as a function of pH at various urea concentrations. We have ignored the data at 4.5 M 
concentration of urea, because there appears to be no significant change in [(f\222 with pH as 
measured by CD spectroscopy. At 3.0 M the value of [(f\222 at the highest pH is considerably 
less than at all other values of urea concentration, which does not seem to be reflected in 



the quantitative analysis presented in Table I of [33J]. Examination of all other curves clearly 
shows evidence for three state behavior. If we calculate Q c for the transition U # N, ignoring 
the equilibrium intermediate I, we find that the overall cooperativity is highly compromised. 
For the purposes of measuring the cooperativity of folding we have focused separately on 
the transitions U # I and I # N and calculated fi c (U # I) and £l c (l # N) using a 
simple generalization of Eq. (|l|). An estimate of <j p h can be easily made by associating 
the midpoint of the transition U # I with the collapse process and the midpoint of the 
transition I ^ N with the folding transition. The values of Q c and a v u for acid induced 
refolding are given in Table III. It is clear that fl c is largest for the smallest values of <j p h, 
which is consistent with the results for LMSC and LM and the two state folders. In Table 
III we also present the cooperativity index f2 c (U ^ I) for apoMb. It is clear from Table III 
that fi c (U # I) is considerably less than Q c (l # N) at all values of urea concentration. The 
very small values of fi c (U # I) suggest that the initial compaction results in the formation 
of an equilibrium intermediate, whose structural characteristics are closer to U than to N. 
Further evidence for our assertion that U # I transition in apoMb is less cooperative (as 



measured by f2 c ) comes from the recent work of Luo et al. These authors measured 

the unfolding thermodynamics from the pH 4 folding intermediate of apoMb, in which A, 
G, and H helices are ordered. By using mutations in the A and G helices Luo et al. argue 
that the unfolding from I is well described by a two state thermodynamics. Although the 
authors did not explicitly provide AGh 2 o values (needed in Eq. (|14|)) we scanned the data 
in Fig. (3) of |37| and calculated fl c values. The resulting fl c values for U # I transition is 
given in Table III. We find that these values are remarkably consistent with those obtained 
from the pH induced unfolding transition U # I reported in . Thus, it appears that the 



assembly of the native state from the equilibrium intermediate I is more cooperative than 
the formation of I itself from the unfolded state. 

Doniach and coworkers have argued using radius of gyration data that folding of lysozyme 
at pH=2.9 is better described by a three state analysis rather than by an all-or-none tran- 
sition [P5] . We have used the theoretical methods developed here to compute a D and Q c . 
These are displayed in Table III. The value of <jd is smaller than for apoMb. Unlike the 
results for apoMb the data for denaturant induced unfolding of lysozyme does not exhibit 
a clear evidence of a third equilibrium intermediate at pH=2.9. This observation together 
with relativity small value of o makes us suggest that at pH=2.9 lysozyme is, at best, a 
borderline three state folder. 

In Fig. (8) we plot Q c as a function of a using the data given in Tables I, II, and III. 
This figure shows a good correlation between two quantities. If the transition is strictly an 
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all-or-none, then from a mathematical representation of the resulting sigmoid curve one can 
show that the correlation between Q c and a should follow independent of the value of C 
used in Eq. ([U]) to get T e (and hence a). Thus, the good correlation between Q c and a for 
two state folders is not surprising. Fig. (8) also shows a good correlation for acid induced 
unfolding of apoMb, which is a three state folder (see Table III). For this protein, Tg and Tf 
are calculated from the midpoints of U # I and I f^N transitions, respectively. Thus, we 
feel that the degree of cooperativity in protein folding can be understood in terms of a with 
Tq being empirically defined by Eq. ([16]). The physically motivated choice for C (< 0.1) is 
further supported by the relationship between folding times and a for several proteins (D.K. 
and D.T., unpublished). 

It can be seen from Table II that there is a good correlation between fl c and AGh 2 o, 
which is the free energy difference between the unfolded state and the native state. The 
value of AGh 2 o is normally obtained by extrapolating to zero denaturant concentration 
which involves errors. More importantly, this analysis of correlation between Q c and AGh 2 o 
cannot be carried out easily in the case of thermal unfolding. As a result it appears that 
it is more convenient to analyze cooperativity data in terms of a, which can be estimated 
using our empirical criterion for both thermal and chemical unfolding transitions. 

By comparing the range of a values for two and three state folders (see Tables I, II, III) 
we can conclude that proteins with a < 0.25 should exhibit an all-or-none transition. If 
a > 0.25, then one should anticipate a three (or more) state folding transition. The value of 
a ~ 0.25, at which crossover from two state behavior to a three state transition takes place, 
should be viewed as a rough estimate. 



IV. CONCLUSION 

One of the most important characteristics of protein folding is that the equilibrium 
folding transition is cooperative ||. The minimal models (lattice and off-lattice) without 
side chains capture this cooperativity only in a qualitative manner [ffl-0]. If a quantitative 
description of this process is required, then, as demonstrated here and elsewhere ||, a coarse 
grained description of side chains becomes necessary. Inclusion of side chains into a lattice 
model of even a small sized chains gives enhanced cooperativity. In these more realistic 
models cooperativity is due to the possibility of denser packing ]16[ of the interior of the 
folded protein. It is clear from the results for N = 27 without side chains that the degree 
of cooperativity approaching that seen in two state folding of real proteins can be easily 
achieved by including side chains. 

We have shown that there is a correlation between the measure of cooperativity Q c (see 
Eq. (|Tl~D) and a = (Tg—Tf)/Tg. Sequences with small a fold in a two state manner with large 
values of Q c . This is the case for lattice models as well as real proteins. The estimates of a for 
two state folders were made using the parameterization of experimental data, measuring the 
fraction of the native state. Since there are two transitions involved it would be desirable to 
obtain Tq using entirely different probe (for example, by measuring the radius of gyration). 
Such measurements would be required to firmly establish the relationship between Q c and 
a. In the absence of such measurements our estimate of a for two state folders based on Eq. 



(|16|) is expected to be only approximate. 
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Proteins, which require a three state model to describe their folding transition, are ex- 
pected to have relatively large a. Based on this classification, we expect that for these 
proteins the actual cooperative transition should occur from an equilibrium intermediate to 
the native state as the temperature or denaturant concentration is altered. Such a coopera- 
tive transition has in fact been recently reported in the folding of apobovine a— lactalbumin 
f38fl . The small values of Q c for the overall U N transition (presumably, resulting from 
large values of a) for lysozyme and apoMb is a reflection of a poor cooperativity of the 
transition from a random coil to the native state due to the presence of the equilibrium 
intermediates. This is consistent with a small values of fl c for transitions U # I, while 
transitions I ^ N are highly cooperative. 

The detailed study of the transition to the native state from unfolded conformations 
in lattice models with side chains suggests a clear mechanism for folding of small proteins 
that exhibit all-or-none transitions. At a relatively high temperature the backbone becomes 
"frozen" adopting a native conformation, but side chains are still mobile. Only at lower 
temperature the side chains become predominantly fixed. The temperature interval, over 
which this happens, is small for sequences with small a, so that these sequences are expected 
to fold kinetically and thermodynamically in a two state manner. This implies that for 
these sequences pathways should exist, in which the acquisition of the native conformation 
(backbone and side chains) and the collapse process should be almost simultaneous. Recent 
molecular dynamics simulations on small peptides also support this notion |[39j| . It is clear 
that the enhanced degree of cooperativity in LMSC over LM is due to packing of side 
chains. A corollary of this result is that denaturation in proteins is basically caused by the 
disruption of the densely packed side chains as was suggested some years ago by Shakhnovich 
and Finkelstein [4(] . 

Sequences with moderate or large a have large difference between the temperatures at 
which the backbone and side chains adopt native-like conformations. As a result these 
sequences get trapped in an intermediate, which inevitably slows down the folding kinetics. 
Explicit kinetic simulations of models with side chains confirm these assertions (D.K. and 
D.T., unpublished). 
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FIGURES 



Fig. 1. Native conformations and energy spectra of sequences A (upper panel) and B 
(lower panel). The sequences for A and B are also displayed. The first letter corresponds 
to residue labeled 1. The energy spectra for the two sequences is given on the right. The 
spectra for the two sequences are divided into two columns. The left column gives the low 
energy levels for conformations, in which the backbone (and hence most likely the entire 
structure) is in a non-native state. The right column corresponds to conformations with the 
backbone in the native state. Sequence A is highly optimized with Z = —17.1 and 5 = 0.50, 
whereas sequence B, while having approximately the same Z score factor of —17.9, is not 
well optimized. This is because one can design a sequence with the same composition as 
sequence B has but with smaller Z. Consequently, the value of 5 = 0.74 for sequence B is 
larger, which indicates a higher degree of heterogeneity of this sequence. Sequence A has 10 
hydrophobic residues, while B has 8. The figure has been generated using RasMol program 
(R. Sayle, 1995). 

Fig. 2. Various thermodynamic functions as a function of temperature for sequence A. 
Fig. (2a) shows < x > (solid line) and < Xbb > (dashed line) (see Eqs. (0) and (|]), 
respectively). The fluctuations A% (solid line) and A\bb (dashed line) are shown in Fig. 
(2b). The temperature dependence of the specific heat C v (solid line) and the derivative of 
the radius of gyration with respect to temperature d < R g > /dT (dashed line) are displayed 
in Fig. (2c). The value of temperature corresponding to the peak in C v is associated with 
T e , which clearly coincides with the peak in d < R g > /dT. This demonstrates that T e is 
the usual collapse transition temperature. The probability of being in the NBA Pnba as 
a function of temperature is given in Fig. (2d) with the horizontal line marking the level 
of 0.5. The temperature, at which Pnba is 0.5, is 0.25, which is almost the same value, at 
which A% has a maximum. Thus, two methods for computing Tf are equivalent. 
Fig. 3. Thermally weighted distribution of states (x, E) calculated for sequence A (upper 
panel) and sequence B (lower panel). The most populated states (x, E) are shown in red. 
This figure clearly shows that, at the folding transition temperature Ty, sequence A samples 
two distinct thermodynamic states - native basin of attraction and the random coil (un- 
folded) states. Sequence B, on the other hand, samples at Tf at least three distinct states 
with different x an d E. We expect the thermodynamic transition for sequence A to be two 
state like, whereas sequence B is better described by a three state analysis. 
Fig. 4 Dependence of cooperativity Q c on <jt for 3D 15— mer sequences without side chains 
(full circles) and with side chains (empty circles). Solid line is a guide to eye. 
Fig. 5 Dependence of cooperativity index Q c on gt for 2D 15— mer sequences with side 
chains and different interaction matrices: random site model (diamonds), MJ shifted poten- 
tials (full circles), MJ potentials (empty circles), KGS shifted potentials (triangles), KGS 
potentials (squares). Shifted MJ or KGS potentials are obtained by setting < B >= 0, 
where < Bq > is the interaction energy averaged over all pairs of residues. The values of ot 
and fl c are obtained by exact enumeration. Solid line is a guide to eye. 
Fig. 6 Dependence of cooperativity index Q c on for 3D 27— mer sequences without side 
chains. The solid line is a guide to eye. 
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Fig. 7 Dependence of the collapse temperature T' e estimated using Eq. fll6|) on actual col- 
lapse temperature Tq determined from the peak of specific heat C v for 27-mer LM sequences. 
The solid line is a guide to eye. 

Fig. 8 Variation of the cooperativity index Q c on a for CI2 (thermal unfolding - full circles, 
chemical unfolding - empty circle), tendamistat (square), and RNase A (triangle), 10 FNIII 
(star), lysozyme (diamond), and apoMb (crosses). We did not show the results for the 
barnase fragment and 9 FNIII, because fl c is too small for them. This figure shows that 
cooperativity can be parameterized by a, which can be experimentally determined. 



20 



TABLES 



TABLE I. Temperature Induced Unfolding: Two State Folders 







experimental values 




theoretical estimates 




Protein 


P H 


T f °C 


T e ° C c 


a 




CI2 a 


2.2 


41.5 


55.5 


0.25 


5.9 




2.5 


47.1 


59.8 


0.21 


9.1 




2.8 


55.1 


67.5 


0.18 


13.1 




3.2 


63.9 


75.2 


0.15 


21.2 




3.5 


70.9 


81.5 


0.13 


29.3 


RNase A b 


8.0 


60.3 


65.2 


0.08 


99.4 



a reported in 
b reported in 
c estimated using Eq. ([16]) 
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TABLE II. Chemically Induced Unfolding: Two State Folders 





experimental values 




theoretical estimates 


Protein 


AGff 2 o, kcal mol -1 


a 




Tendamistat a 


8.13 


0.17 


13.7 


CI2 6 


7.03 


0.20 


9.9 


10 FNIII c 


6.1 


0.22 


7.5 


9 FNIII c 


1.2 


0.60 


0.28 


Barnase fragment^ 


1.1 


0.62 


0.25 



a GdnHCl induced denaturation reported in |H 
b GdnHCl induced denaturation reported in | 
c GdnHCl induced denaturation reported in [^2| 
d Trifluoroethanol induced denaturation reported in |33 
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TABLE III. Three State Folders 



Protein 


P H 


C, [Urea] 


a 


o c (u # i) 


n c (i # n) 


Lysozyme a 


2.9 




0.21 




9.1 


ApoMb 6 




0.0 


0.41 


3.4 


10.5 






1.0 


0.34 


4.3 


19.3 






1.5 


0.32 


3.7 


24.7 






2.0 


0.32 


2.7 


40.2 


ApoMb c 


4.2 






1.2 d 






4.2 






5.0 e 





a Urea induced denaturation at pH=2.9 and T = 20°C monitored by R g reported in [3^| 
b pH induced denaturation at T = 0°C monitored by CD signal reported in |34} 
c Urea induced transition at pH=4.2 and T = 4°C monitored by fluorescence reported in 
37 

d 20 mM Na 2 SO A 
6 50 mM NaClO A 
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